

float cor(float mx, float my, float sdx, float sdy, __global float *cx, __global float *cy,  int num)
{
	int i;
	float sum = 0, cov;
	if (num == 0)
		return 0;
	else
	{
		for (i=0; i < num; i++)
		{
			float difx = (cx[i]- mx);
			float dify = (cy[i] - my);
			sum = difx*dify + sum;
		}
	}
	cov = (sum / num);
	return cov/(sdx*sdy);
}

float mean(__global float *m, int num)
{
	int i;
	float sum = 0;
	if (num == 0)
		return 0;
	else
	{
		for (i=0; i < num; i++)
			sum = m[i] + sum;
	}
	return (sum / num);		
}

float sd(float means, __global float *s, int num)
{
	int i;
	float sum=0;
	if (num == 0)
		return 0;
	else
	{
		for (i=0; i < num; i++)
			sum = (means - s[i]) * (means - s[i]) + sum;
	}
	return sqrt((sum / (num - 1)));
}

float pearson(__global float *xx, __global float *yy, int num)
{
	float meanx = 0, meany = 0, stddevx = 0, stddevy = 0;

	meanx = mean(xx, num);
	meany = mean(yy, num);
	stddevx = sd(meanx, xx, num);
	stddevy = sd(meany, yy, num);
	return cor(meanx, meany, stddevx, stddevy, xx, yy, num);
}

__kernel void MutualInfo(__global float *info, __global float *arr, int M, int NUM_ROWS)
{

	const int idx = get_global_id(0);

	__global float *pRowx, *pRowy;
	pRowx = &arr[idx * M];

	if (idx < NUM_ROWS)
	{
		int jdx;
		for (jdx = idx; jdx < NUM_ROWS; jdx++)
		{
	
			pRowy = &arr[jdx * M];
			info[idx * NUM_ROWS + jdx] = pearson(pRowx, pRowy, M );
			info[jdx * NUM_ROWS + idx] = info[idx * NUM_ROWS + jdx];
		}
	}
}


